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We study the surface roughness of prototype models displaying self-organized criticality (SOC) 
and their noncritical variants in one dimension. For SOC systems, we find that two seemingly 
equivalent definitions of surface roughness yields different asymptotic scaling exponents. Using 
approximate analytical arguments and extensive numerical studies we conclude that this ambiguity 
is due to the special scaling properties of the nonlinear steady state surface. We also find that there 
is no such ambiguity for non-SOC models, although there may be intermediate crossovers to different 
roughness values. Such crossovers need to be distinguished from the true asymptotic behaviour, as 
in the case of a noncritical disordered sandpile model studied in [10]. 



Since the original proposal by Bak, Tang and Wisen- 
feld [1] there has been a large body of work directed to- 
wards understanding ubiquity of scale invariance in exter- 
nally driven open dissipative systems using the concept 
of self-organized criticality (SOC). The sandpile is a pro- 
totype model system which has been extensively used as 
a paradigm of SOC [2,3]. The principal aim has been to 
elucidate how slowly driven dissipative systems with fast 
relaxation mechanisms display long-tailed distributions 
of activity sizes. Recently, there have been resurge of 
interest in sandpile models in order to understand SOC 
in connection with better understood scale invariance in 
other noncquilibrium phenomena such as absorbing state 
phase transitions [4] and driven interfaces [5]. 

In [6], Krug, Socolar and Grinstein (KSG) studied the 
surface fluctuations in a prototype model of SOC, the 
limited local sandpile (LLS) [3] and its variations in one 
spatial dimension and concluded that the interfacial fluc- 
tuations, although nontrivial, are evidently unconnected 
to the criticality of the system. In fact, they argued 
that the evolution of height fluctuations h(x,t) can be 
described by an extension of the Kardar-Parisi-Zhang 
(KPZ) [7,8] equation for an anchored interface 



dth = Ddlh + cdjh + X{d x hf + rj(x, t) , 



(1) 



where T)(x, t) is a Gaussian white noise. The linear 
term cd x h is more relevant than the diffusion term Dd^h 
and the KPZ nonlinearity X(d x h) 2 (in RG sense) and, 
since the interface is anchored, can not be eliminated 
by a Galilean shift. This term is responsible for trans- 
porting fluctuations up the pile thus relating the spa- 
tial roughness of the anchored interface to the dynam- 
ical roughening of a moving KPZ interface in d = 1, 
otLLS — I^KPZiid) = 1/3- This argument can be easily 
extended to higher dimensions. 
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FIG. 1. (a) Dynamics of the sandpile: An initial surface 
(thick line) on becoming unstable relaxes by transferring the 
unstable cluster (shaded blocks) successively to the left till 
a stable configuration is reached, (b) Log-log plot of the 
avalanche size distribution P(m) for the LLS (solid symbols) 
and the ILLS (open symbols) for two different system sizes. 
P(m) for the ILLS falls off with a characteristic length scale 
independent of L. 

In this paper we study interface roughness in LLS and 
related models and find that effect of the criticality of the 
system is somewhat subtle and is due to the nontrivial 
nature of the steady state surface. The value of surface 
roughness depends on whether it is measured about this 
nontrivial steady state profile or about the mean instan- 
taneous surface, the latter being the correct choice for 
moving interfaces. The assertion made by KSG is valid 
only if h in (1) is defined as the fluctuation around the 
nonlinear steady state surface. We also study modifica- 
tions of the LLS which have nonlinear steady state pro- 
files but which do not display SOC and find that there 
is no such ambiguity in the asymptotic roughness expo- 
nent. Thus, we conjecture that the roughness exponent 
is uniquely defined for systems not displaying SOC even 
though they may possess nontrivial steady state surfaces 
[9]. For SOC systems the ambiguity exists and one needs 
to define the roughness appropriately. Our results have 
important bearings on studies such as the recent one by 
Barker and Mehta [10] who observed anomalously large 
roughness exponent in a sandpile model with structural 
disorder. 



The limited local sandpile model in d = 1 is defined 
as follows (see Fig. la): on a one dimensional lattice of 
sites i — 0, 1, • • • , L we define an integer height variable 
hi. One grain of sand added at a randomly chosen site 
increases the height at that site by one: hi — ► hi + 1. 
The configuration is stable if the local slopes satisfy 
Zi = h. i+ i — hi < z c for all i, where we chose z c = 2 
(our results are essentially unchanged for z c > 2). An 
instability occurs when by addition of grains at site i + 1, 
the local slope exceeds threshold Zi > z c ; in this case 
z c grains are transferred from column i + 1 to column i. 
Subsequently, columns i and i + 2 may become unstable, 
setting off further topplings leading to an avalanche. The 
grains involved in an avalanche leave the system if they 
reach site i = (ho(t) = 0,Vi). A new grain is added af- 
ter the system has attained a stable configuration. The 
final configuration reached is independent of the order in 
which the sites are updated in case more than one site be- 
come unstable. We count one unit of time (Monte Carlo 
step) for every L grains added. 

The standard way of quantifying avalanches in these 
systems is to count the number m of grains that leave 
the pile after each deposition [3]. The criticality of the 
model is reflected in a broad probability distribution 
of avalanche sizes, P(m) (Fig. lb), which has no other 
length scale except the system size L. 

In [6], KSG also introduced and studied the inertial 
limited local sandpile (ILLS) which mimics the effect of 
inertia of the falling grains in a real sand pile. The insta- 
bility condition setting off an avalanche is same as that 
in the LLS but the condition of stopping is changed. A 
cluster of grains when first destabilized is assigned a mo- 
mentum p = and each time the front of the cluster 
reaches an unstable site (zi > z c ) it gathers momentum, 
p — ► p + 1. If it comes across a stable site, p decreases by 
an amount r: p — > p — r. The cluster continues moving 
down as long as p > and leaves the pile if it reaches 
the site i = 0. Clearly, r = oo corresponds to the LLS. 
It was noted by KSG that ILLS is not critical for any 
finite value of r which is reflected in the corresponding 
avalanche size distribution (Fig. lb). 

To study the interfacial fluctuations in these models, 
we start with an initially flat pile (hi(0) = 0, Vi) and add 
grains till the pile reaches a steady state with the mean 
surface making a critical angle ip (taaipL-^oo = 3/2) with 
the horizontal. The closed boundary condition at i = L 
ensures that the pile has only one surface with the crit- 
ical slope [9]. The width of the interface, in the steady 
state, can be measured in two ways 

W?{L) = ^Y,(Mt)-(hi)Y) and, 

^2 2 (£) = ^£«)-s(*H] 2 >, (2) 

i=l 



where s(t) = 2[£ti h t (t)]/L(L + 1) is mean slope of the 
interface about which sum of instantaneous height fluctu- 
ations vanish [h>i{t) — s{t)-i\ =0. In both definitions 
(2) above the ensemble average (• • •} is identical to the 
time average in the steady state. 
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FIG. 2. Width vs system size for LLS: W\ (circles) and W2 
(squares) are obtained using (2) while Ws (triangles) is com- 
puted from (5) below. W c (diamonds) ~ L 018 asymptotically. 
Inset: W\ (circles) and W2 (squares) and Ws (triangles) for 
the ILLS with r = 20. 

For large enough system sizes the scaling hypothesis 
is expected to be valid and the roughness exponent a is 
defined through 

W h2 (L)~L a , L^oo. (3) 

In Fig. 2, we plot the two widths W\ and W2 for a 
range of system sizes L = 2 n with 6 < n < 14 (that 
is 64 < L < 16384). It is seen that the two widths have 
different asymptotic behaviours and give two different 
values of a. The roughness a ~ 0.33 computed from W± 
(and W2 up to L < L c ) is in accordance with the pre- 
dictions of Eq.(l). The asymptotic roughness computed 
from W2, beyond the crossover system size L c ~ 1024, 
a = 0.65 ± 0.02 is larger than that for an interface gener- 
ated by a simple random walk. This is surprising since, 
prima facie, the definitions in (2) are expected to be 
equivalent as far as asymptotic scaling is concerned ( 
[11]). In fact, for the ILLS, we have computed W\ and 
W2 for r = 20 (Fig. 2, inset) and it is seen that they 
indeed have the same asymptotic behaviour and hence, 
a unique value of a ~ 0.33. 

We now show that the crossover in W 2 m the LLS can 
be traced to the fact that the time averaged steady state 
profile is not linear, i.e., (hi(t)) ^ s i and this makes the 
two definitions nonequivalent. Although, this nonlinear- 
ity of (hi(t)) is a consequence of the boundary effects and 
is present for the ILLS as well, we demonstrate that the 
presence of SOC in LLS results in special scaling prop- 
erties which is responsible for the observed differences 
between asymptotic scalings of W\ and W2. It can be 
shown that W\ and W2 are related by 



W 2 (L) + W 2 (L) = W 2 (L) + W§(L) . 



(4) 



where Ws(L) is the root-mean-square (rms) wandering 
of the steady state interface profile (hi), 



Wl(L) 



S ■ l\ 



(5) 



i=0 



with s = (s(t)) being the mean slope of the steady state 
profile satisfying ^2 i=0 [(hi) — s ■ i] = 0. The second term 
on the left of (4), W 2 {L) = ([s(t)-s] 2 )(L+l)(2L+l)/6, is 
the contribution to interface width due to instantaneous 
slope fluctuations and is the analogue of the center of 
mass fluctuations in the case of moving interfaces [8] . In 
fact, in the case of moving interfaces W c dominates for 
large length scales and hence W2 is taken as the width. 
As it turns out, for the anchored interfaces we deal with 
here, W c is negligible compared to the other terms in (4) 
as L — > 00. From Fig. 2, it is evident that the crossover 
in Wi{L) for the LLS is due to the contribution from Ws- 
In the following we show that, for the LLS, the special 
form of the nonlinear steady state profile (hi(t)), which 
results in Ws ~ L 66 , is a consequence of the singular 
diffusion associated with the SOC state. 

In [12], using a simple model, Carlson et. al. demon- 
strated that the SOC state of the system is associated 
with the vanishing of the density of troughs concommit- 
tant with the divergence of the corresponding diffusion 
coefficient. The troughs are defined as the sites for which 
Zi < so that, in the LLS, an avalanche necessarily stops 
at a trough (or leave the system at i = 0). On a coarse 
grained level we define a set of densities {p n (x,t):n = 
—00, • • • , — 1, 0, 1, 2} where p n (x, t) denotes the local den- 
sity of sites with zi = n. It follows that the coarse 
grained local slope z(x,t), which is locally conserved by 
the dynamics, may be expressed in terms of the p„'s as 
z(x,t) — ^2 n —_ 00 np n (x,t). The open boundary condi- 
tion at i = implies /i(0, t) = and the closed boundary 
condition at i = L is modeled by setting z(L, t) = 0. The 
trough density p(x,t) = J2n=-oo Pn( x , t) is not strictly 
conserved. In the '01' model considered in [12] the slope 
has only two values Zi = (trough), 1 and hence both z as 
well as p are strictly conserved and are related simply as 
z = 1 — p. Although p is not strictly conserved for LLS 
and ILLS we still approximate its dynamical evolution 
by the continuity equation, d t p(x,t) + d x J(p(x,t)) = 0. 
Phenomenologically, the trough current J consists of 
three parts: (i) current due to addition of grains Jo, (ii) 
avalanche current J a, and (iii) a microscopic noise term 
rj(x,t). In analogy with driven diffusive systems [6] the 
phcnomcnological form of the avalanche current is writ- 
ten as Ja{p) — ap + bp 2 + ■ ■ ■ — D(p)d x p. As in the 01 
model, the critical state of the system (L — > 00, r = 00) 
is associated with p, d x p —> and hence in order to 
balance the finite input flux Jo, the diffusion constant 
D(p — > 0) has to diverge appropriately [12]. Without 



loss of generality the leading divergence in D is taken to 
be a simple pole and the evolution of p(x,t) is thus [13], 



dp d 
dt dx 



D(p) 



A(p) 



(6) 



Here A(p) is taken to be a smoothly varying function in 
the relevant interval < p < 1. 

In order to compute the steady state density pro- 
file p(x) = p(x, t — > 00) of troughs, wc note that, as 
L — > 00, the system is arbitrarily close to criticality 
(p — ► 0), and hence A(p) may be set to a constant 
A = A(0). One can then readily integrate the current 
equation D(p)dp/dx — — J , subjected to the boundary 
condition p(L) — p^, to obtain 



p(x)^ PL [I + 7 (L-x)Y 



(7) 



where = 1/(0- 1) and 7 = -J /(M ). The aver- 
age density of the troughs p = L^ 1 p(x)dx thus scales 
with system size asymptotically as p ~ pLi^L)^ 6 / (1 — 6). 
From Fig. 3 (inset), we obtain 6 ~ 0.33 and thus <f> ~ 4 
for the LLS, which are the same as the corresponding 
values obtained in [12] by direct measurement of D(p) in 
a closed system. 

In the following we find an approximate numerical re- 
lation between p(x) to z(x) which would enable us to 
find the steady state interface profile ho(x) — J Q X z(x)dx 
form (7). In Fig. 3 (inset), we plot the average den- 
sities ~p n = L^ 1 J a dx(p n (x,t)) in the steady state as 
functions of system size L. We note that, as L in- 
creases, while J) 1 ~ ~p 2 ~ 1/2 remain finite, densi- 
ties of all the trough sites vanish algebraically with L: 
PojP-i ~ i _0,33 ,P-2 ~ £~ ' 6 , and p n 's with n < — 3 are 
negligible. Hence, in the limit of large L, we can approx- 
imate the total trough density as p ~ p + ~p~- 1 • From 
the normalization condition X)n=-oo Pn( x , t) = 1 it fol- 
lows that pi + p2 = 1 — p and we find numerically that 
Pj — p 2 ~ 0.46p. Thus, for large systems, we may write 
z (x) ~ 2p 2 + pi — p-i ~ 3/2 — np(x), where numerically 
k ~ 2 [14]. Thus, ho(x) is given approximately by 



h (x) 



Kp^l+jL) 1 - 

J(l-0) 



1 + 7(L - x) 
1+jL 



(8) 



In Fig. 3, using the approximate expression for p(x) from 
(7) with (f> = 4, we compare ho (x) with that obtained nu- 
merically and notice that the agreement is rather good 
given the nature of approximations involved [15]. 
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FIG. 3. Deviation of the steady state height profile from 
straight line. For the LLS, the simulations data (solid curve) 
is well approximated by the analytic form (8) (dashed curve) 
for L = 2048 with optimally chosen values of parameters pl 
and 7. The dotted curve is for the ILLS with r = 20. Inset: 
Average density ~p~ n of sites with slope Zi = n in the steady 
state as a function of system size L. 

The L dependence of Ws of the steady state profile 
can be computed from (8) and turns out to be 



W S (L) -L 1 - 



(9) 



For the LLS with 9 ~ 0.33, W S (L) ~ L om dominates 
Wi ~ L 1 / 3, and W c ~ i 018 in (4) and hence, for large 
i, W2(L) ~ i - 66 in very good agreement with our 
numerical results (Fig. 2). It is interesting to note that 
W\ and W2 would have different asymptotic behaviours if 
9 < 2/3, i.e., if > 5/2. For the '01' model studied in [12] 
it is shown that <j> = 3 exactly and thus Wi{L) ~ L 1 ! 2 for 
the corresponding interface [16]. Our conjecture that W\ 
and W 2 would have different asymptotic behaviour for 
SOC systems would be invalid if one can devise a model 
with 4> < 5/2. 

Next, we present a general argument as to why, for 
noncritical models such as the ILLS, wc do not expect 
any ambiguity in roughness exponent. The essential dif- 
ference between a critical and a noncritical model is the 
presence of an additional length scale (apart from system 
size and the microscopic cutoff) which shows up in the 
avalanche size distribution (e.g. Fig. 1 for ILLS) as well 
as governs the decay of the boundary effects into the bulk. 
In the ILLS this length scale is related to the parameter 
r. To see this we first note that as a cluster moves down 
its momentum p makes a random walk [6] . If the mean 
density of troughs p < 1/r then most avalanches leave 
the system resulting in net drainage and if p > 1/r then 
most avalanches stop on the pile leading to net growth 
of the pile. Thus, in the steady state one has a finite 
density of troughs p — 1/r and thus the mean spacing 
between the troughs 1/p = r appears as the additional 
length scale. As already discussed above, the current 
J(p) now has a nonzero systematic part ap + bp 2 . . ., in 
addition to the finite diffusion term. Such a systematic 
current, e.g., a term such as ap, results in the effects of 
the boundaries decaying exponentially inside the bulk. 



This is in contrast to the long ranged power law decay in 
the LLS which is reflected in the forms of p(x) and ho(x). 
Although the steady state surface is nonlinear (Fig. 3), 
its rms wandering Ws is bounded and does not alter the 
asymptotic scaling of W 2 ■ 

Lastly, we briefly point out possible pitfalls of using 
W2 naively without properly accounting for the under- 
lying steady state profile, even in systems which do not 
show SOC. Recently, Barker and Mehta [10] studied a 
disordered version of the LLS where disorder was intro- 
duced by allowing the added grains to have an aspect 
ratio different from unity: grains are now rectangles and 
are deposited on the sandpile with fixed probabilities of 
landing on their larger or smaller edges. Thus, the height 
hi of column i, no longer an integer, is the sum of the ver- 
tical edges of all the grains in that column. In addition to 
the threshold dynamics of LLS, dynamical reorientation 
of cluster of grains were allowed. They found that the 
larger sandpiles cease to display SOC which is reflected in 
the emergence of a preferred size of the large avalanches 
in the associated drop number distributions. Their nu- 
merical studies of W2 showed that while for very small 
systems [L < 100) the roughness exponent a ~ 0.34, it 
seems to crossover to a much larger value a ~ 0.72 for 
larger system sizes, 100 < L < 400 (Fig. 4). Our pre- 
liminary numerical studies, using the same parameters 
as in [10], of yet larger systems (500 < L < 2048) show 
that in fact the crossover seen in [10] is evidently tran- 
sient and the asymptotic roughness exponent goes back 
toa~ 0.33 (see Fig. 4). 
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FIG. 4. Interface width of the disordered sandpile as a 
function of system size. The squares represent W2 and circles 
represent W = \JW% — W§. The parameters chosen are the 
same as in [10]. The straight lines have slope 0.33 (dotted) 
and 0.72 (dashed). Inset: The steady state interface profile, 
after subtracting the mean linear profile. 



[10] is less notice- 
instead of W2 



In fact, the crossover reported in 
able if one looks at W — ^JW% 
(Fig. 4) [17]. This is in accordance with our conjecture 
since the disordered system is not critical, although the 
steady state profile is nonlinear (Fig. 4, inset), it docs not 
affect the asymptotic roughness measured by W 2 - The 



transient crossover seen in W2 here is similar to what is 
seen for the ILLS in Fig. 2, inset. 

In summary, we have studied the surface roughness of 
a prototype model of SOC and its modifications in one 
dimension. We find that one needs to be careful in defin- 
ing quantities such as the interface width since special 
form of the steady state shape of the surface in systems 
with SOC can result in different asymptotic behaviours 
of otherwise equivalent definitions. Although there is no 
such ambiguity for noncritical models, still there may be 
crossovers at intermediate length scales which should not 
be taken as the true asymptotic behaviour. 
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